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ABSTRACT 

We are developing a 3D radiation hydrodynamics code to simulate the interaction of convec- 
tion and pulsation in classical variable stars. One key goal is the ability to carry these simulations 
to full amplitude in order to compare them with observed light and velocity curves. Previous 2D 
calculations were prevented from doing this because of drift in the radial coordinate system, due 
to the algorithm defining radial movement of the coordinate system during the pulsation cycle. 
We remove this difficulty by defining our coordinate system flow algorithm to require that the 
mass in a spherical shell remain constant throughout the pulsation cycle. We perform adiabatic 
test calculations to show that large amplitude solutions repeat over more than 150 pulsation 
periods. We also verify that the computational method conserves the peak kinetic energy per 
period, as must be true for adiabatic pulsation models. 



Subject headings: convection — hydrodynamics 
variables: general — stars; variables; RR Lyrae 

1. INTRODUCTION 

Early nonlinear calculat i ons o f stellar pulsa- 
tion, as outlined bv lChristvl (|l964l) . used a ID La- 
grangian framework and had considerable success 
producing full amplitude RR Lyrae models that 
resembled the basic observations. However, these 
calculations used purely radiative envelopes and 
failed to identify a red edge to the RR Lyrae insta- 
bility strip. This lead to the hypothesis that con- 
vection in the ionization zones quenched p ulsation 
(|Baker & Kippenhahn"l965l;'ChTistyll966!). To ex- 
plore this, Tugglc & Ibcn (197^ used a ID, linear, 
non-adiabatic code with a time-independent mix- 
ing length theory for the convective flux. They 
found that time-independent convection only re- 
duced the growth rate of pulsation but did not 
produce pulsational stability. 
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stars: 



Several formalisms and ID codes were devel- 
oped that included time-dependent convection 
by introducing an additional differential equation 
to calculate the evolution of the convective flux 
with ti me based on the s ta ndard mixin g length 



theor y (iCox et al 



1966allbl: lUnnd Il967l: iGoueh 



Il977l) . IStelhngwerd (Il982allbl . [l984allbl ig) developed 
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a time-dependent treatment of convection for ID 
Lagrangian models which follows the time evolu- 
tion of averaged convective velocities and includes 
a treatment of overshooting and eddy viscosity 
to account for small length scale kinetic energy 
dissipation. 

There are several numerical difficulties associ- 
ated with modeling radial pulsation. First there 
are very steep gradients in the ionization zones 
and adequately resolving thes e gradients is impo r- 
tant for a ccurate modeling. Gehmeyr (Il992al lbl. 
19931) and iDorfi fc Feuchtingeil (|l99ll) have both 
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included adaptive grids too better resolve the 
steep gradients in the ionization zones as they 
sweep through the envelope during pulsation. Us- 
ing a version of Stellingwerf 's time-depe ndent con- 
vective model with his adaptive scheme, iGehmevi 
was able to produce a red edge at roughly the ob- 
served effective temperature. He notes that the 
effective temperature of the red edge is dependent 
on the parameters used for the convective model, 
and that the predicted temperature of the red 
edge could vary by a few hundred degrees Kelvin 
depending on the values used for the convective 
model parameters. Also, there are differences be- 
tween Gehmeyr's light amplitude-rise time rela- 
tionship and the obs erved relationship in both 
slope and zero point. Feuchtinger fc Dorfi ( 1996[ ) 
used their adaptive code to calculate light and ra- 
dial velocity curves which exhibit shapes typical 
of RR Lyrae stars. A second potential difficulty is 
an accura te representation of t he surface bound- 
ary, which lFbuchtinger fc Dorfi tested by including 
a model atmosphere and found that its inclusion 
did not impact the pulsational characteristics of 
the model. 

iMarconi et al. ( 2003f ) used the ID, Lagrangian, 
hydro dynai nics code de s cribed Ijy iBono fc Stellingwerf 
(|l994l ) and iBono et all (|l997allbh to compute RR 
Lyrae models to compare with the RR Lyrae stars 
observed in M3. In order to fully specify the 



problem lMarconi et al.l needed to choose a mixing- 
length parameter, and adopted both I /Hp = 1.5 
and 2.0, where I is the mixing-length and Hp is 
the pressure scale height. They found that in or- 
der to match the boundaries of RR Lyrae gap in 
M3, they required two different mixing-length pa- 
rameters, one to obtain the observed blue edge 
location {I /Hp w 1.5) and the other to obtain 
the observed red edge location {l/Hp w 2.0). 
In addition the observed visual amplitude as a 
function of B-V displays nonlinear characteristics, 
while theoretical rela tions predict linear relation- 
ships. iMarconi et al.l also mention that a mixing- 
length parameter of 2.0 produces luminosities for 
horizontal-branch models that are brighter than 
what is observed by ss 0.08 ± 0.05 mag. 

Other models for time-dependent conv ection in 
one dimensi on ha. v e bee n proposed by iKuhfuss 
(|l986,) and .Xiond (Il989l). Kuhfuss argued that 
the convective model bv lStellingwerfl (|l982al ) does 
not use the diffusion approximation consistently 



throughout the model. ISmolec fc MoskalikI (|2008[ ) 
developed their application of the KuhfussI convec- 
tive model, which requires eight free parameters, 
and used it to study conv ection in /? Cephei stars 
(|Smolec fc MoskalikI [2OQ7I ). They found that con- 
vection is not important for calculating pulsation 
amplitudes for their models. However, they cau- 
tion that their convective model, while working 
well for classical pulsators, is at the limits of its ap- 
plicability in the /3- Cephei models they are s tudy- 
ing. More recently, lOlivier fc Woo d' ("2005!) have 
also developed a code including the Kuhfus^ con- 
vective model and present some test calculations 
of their program, mentioning that the turbulent 
viscosity parameter shows potential as an impor- 
tant determinant of the pulsation amplitudes. 

The distillation of multi-dimensional convec- 
tive phenomenon to one dimension is always ac- 
companied by extra equations and/or parameters 
to approximate the effects of convective motions 
of mater i al in m ore than one spatial dimension. 
Deupred (|l977al) approached the interaction of 



convection and stellar pulsation in a fundamen- 
tally different way using a two-dimensional hy- 
drodynamic code directl y following the convec- 
tive flow patterns. While lDeupre3 (|l977bll3 . [l980l 
ll985 ) was able to successfully determine the ob- 
served edges of the RR Lyrae instability strip, he 
was unable to compute full amplitude solutions 
because his algorithm for moving the radial coor- 
dinate allowed the radial zoning to drift over time. 
Consequently at later times, the radial zoning did 
not cover the hydrogen ionization zone adequately 
and the calculations were eventually numerically 
unreliable. 



The algorithm iDeupred used for the 
moving radial coordinate used the horizontal av- 
erage of the radial velocities at a particular radius 
as the grid velocity. Recently Bruenn et al. ( 2010[ ) 
have used a similar average radial velocity as the 
grid velocity to follow the core in fall phase in su- 
pernovae simulatio ns. An other form of a radial 
moving grid was bylMundprecht (2001), where the 
radial grid velocity at the surface was set as the 
average of the radial velocities, and the inner grid 
velocity was held constant. The intermediate ra- 
dial grid velocities were then set using a dilatation 
factor. 

Recently IStokll (l2008l) developed an approach 
for mode l pulsatio n somewhere between the 2D 
model of lDeupr"ee and the ID convective models. 
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Stokll used two radial columns to model convec- 
tion. One column represented the sum of all up- 
ward convective flows, and the other column rep- 
resented the sum of all downward flows. While 
not including any mixing length parameter, it does 
contain free parameters related to the physical size 
of the convective cells, and the fraction of the 
surface area of a spherical surface which contains 
downdrafts. These parameters do have a physical 
basis, but in practice it would be difficult to deter- 
mine them as they are probably functions of depth 
and likely depend on a particular star's proper- 
ties. Also, because the model only has two radial 
columns it may miss some of the more subtle fea- 
tures of convection. More recently others have be- 
gun working on directly simulatin g the interaction 



of convection and pulsation in 2 D (jMuthsam et al 
2OIOI : iGastine fc Dintransll2010l ). 



Both iBuchled (|2009l ) and iMarconil (|2009l ) have 
stressed the importance of improvi ng the co nvec- 
tive models used in variable stars. iBuchlerl high- 
lights some of the well known difficulties facing 
time-dependent mixing length noting that it is 
an empirical description, rather than a consis- 
tent physical description. Also, the up to eight 
or more free parameters used in time-dependent 
mixing length approach can not be chosen based 
on physics, but instead must be calib rated by 
comparison with observations. Marconil mentions 
some remaining problems for RR Lyrae models, 
particularly the unsatisfactory match between the- 
oretical light curve morphology and observed light 
curve morphology near the red edge of the RR 
Lyrae instability strip, supporting the suggestion 
that an improved treatment of turbulent convec- 
tion is needed. 

3D convective simulations have had significant 
success i n other areas of stellar astrophysics. For 
example iNordlund et al ] liooa) note many of the 
recent successes in 3D modeling of solar surface 
convection, in particular that 3D models with nu- 
merical resolutions around 200'^ reproduce widths, 
shifts, and shapes of observed photospheric spec- 
tral lines with high accuracy . 3D c onvective simu- 
lations bv lMeakin fc Arnetti (j2007l ) simulated core 
convection in a massive star finding differences in 
2D and 3D convective velocities in both morphol- 
ogy and magnitude. 3D models remove the need 
for many free parameters and the lack of a phys- 
ically consistent description of convection, with 



convection resulting naturally from the conserva- 
tion laws; however, an algorithm to include the 
effects of sub-grid scale turbulence on the larger 
ed dies is sti ll requir ed. Here we build on the ideas 
of lDeupre3 ( 1977a ). with the goal of developing 
a fully 3D calculation in which the radial coordi- 
nate moves in such a manner to allow us to per- 
form full amplitude solutions of RR Lyrae models 
with Large Eddy Simulations (LES) of convective 
energy transport. As a first step we apply our ap- 
proach to purely adiabatic models to verify that 
the method can compute accurate and stable large 
amplitude periodic solutions over many periods. 

2. HORIZONTAL EULERIAN RADIAL 
LAGRANGIAN SCHEME 

Calculation of full amplitude solutions requires 
following the pulsation for many periods. ID 
codes have bee n able to calculate full amplitude 
solutions, while iDeupred 's 2D code had difficulty 
after many periods because of his particular mov- 
ing radial coordinate. This led us to try using 
the internal mass, Af^, as the radial independent 
variable instead of radius, r, and allow r to change 
such that the mass within a shell remains constant. 
This requires introducing a grid velocity, wor, in 
the radial direction that dictates how the coordi- 
nate system radius changes. The intent is that our 
radial grid acts like that of a ID Lagrangian code 
while allowing the usual Eulerian approach in the 
horizontal directions. Note that this approach still 
allows fluid flow across radial zone boundaries. It 
just moves the radial gridding so that it maintains 
the mass in a spherical shell. It does not put any 
constraints on the horizontal flow and does not 
alter the physics of the conservation equations in 
any way. 

The calculations are performed in a limited 
range of the spherical polar coordinates 9 and 0, 
a 3D version of a "pie slice". Periodic bound- 
ary conditions are placed in the horizontal direc- 
tions. The interior boundary is placed at a loca- 
tion deeper than 0.15 of the stellar radius and is 
regarded as rigid. This is a common assumption 
in most ID simulations. Because the horizontal 
zoning would get very narrow (leading, through 
the Courant condition, to undesirably short time 
steps) and because the 3D flow of interest is ex- 
pected to be only in the surface ionization regions. 
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we impose a purely radial region for an arbitrary 
number of radial zones above the interior bound- 
ary. In conjunction with the assumption that non- 
radial motion occurs only near the very low mass 
surface, we assume the gravitational force has its 
spherically symmetric form. 

2.1. Conservation Equations 

We first define a horizontally averaged density 
which allows us to replace an infinitesimal radial 
change, dr, with an infinitesimal internal mass 
change, dM^. The mass of a spherical shell of 
thickness dr is given by dMr = 47rr^(p)dr where 
(p) is the volume averaged density within a spher- 
ical shell. 



' ,k^i,j,k • 



(1) 



Here i, j, and k are the f , 9 and zone indices 
defined at zone centers and increase with each of 
their respective spherical coordinates (e.g., i in- 
creases from the center of the star towards the 
surface). Vi is the volume of a spherical shell at a 
particular radial shell, i, that spans all the j and k 
at that i. Vij^k is the volume of the {i,j,k) grid cell. 
In order to develop and test this approach, we have 
assumed that the system is adiabatic. Introduc- 
ing both dMr and the grid velocity, vor, into the 
3D conservation equations for mass, three compo- 
nents of momentum, and energy produces: 
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The above symbols have their usual meaning [E 
is the specific internal energy) . Finally the system 
is closed with three further equations: 



dr 

9^ = "°'^ 



which is used to update the radius, 
P = (7 - l)pE, 



(7) 



(8) 



a simple gamma-law gas for the equation of state, 
and an equation for solving for the grid velocity 
(see equation [T5)) . Before we present the equa- 
tion for the grid velocity, we present an equivalent 
equation to equation Q by which we solve mass 
conservation. 

2.2. Finite Volume Mass Conservation 

The definition of volume we will use in the equa- 
tion for determining the grid velocity (eq. [T5] ) 
and the average density (eq. [T]) must be con- 
sistent with the mass conservation equation. The 
definition of volume we use in equations ()15p and 
© is 



riJrl/2 rOj + l/2 r<Pj + l/2 



1/2-' 6j- 1/2'' "/"i- 1/2 



r^ sin 9 d<j)d9 dr. (9) 



To make the mass equation consistent with this 
definition of the volume, one integrates equation 
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Fig. 1. — Geometry of a cell in spherical coordi- 
nates. The areas of the six surfaces of the "cube" 
are shown. 

([2]) over the volume, then uses Gauss's theorem 
to convert the volume integral into a surface inte- 
gral, producing the finite volume form of the mass 
conservation equation. 

The mass in a cell changes only from mass flow- 
ing into and out of that cell. This change from one 
time step (n) to the next {n + 1) can be written as 
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where A is the area of a cell face (see figure [T|) and 
the Fs are fluxes defined as 
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To obtain the densities at the interfaces (p"j_i/2) 
straight averages are computed of centered densi- 
ties adjacent to the interface. The two centered 
subscripts have intentionally been omitted in the 
expressions for the fluxes and areas to reduce equa- 
tion length (e.g. i was left off but i + 1/2 was 
kept). Note that the densities in the fluxes are 
at time n and not n -I- 1/2, so our solution algo- 
rithm is straightforwardly explicit, as is true for 



many ID calculations. It is not possible to prop- 
erly time center all terms without introducing a 
more complex implicit or multi-step explicit algo- 
rithm. With this expression for the fluxes, we can 
directly solve equation pUj) for the density at the 
new time step. 

2.3. Radial Grid Velocity 

The final piece needed to complete the descrip- 
tion is the calculation of the grid velocity. For a 
spherical shell to have constant mass, the net flow 
of mass into and out of that spherical shell must 
be zero. Summing up all the fluxes into and out of 
the individual horizontal cells in a spherical shell, 
substituting equation ([TT|) in for the outer radial 
flux (at i -f 1/2) and setting the result equal to 
zero we arrive at the equation. 
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Solving for the outer grid velocity, v^, 
duces an equation for calculating the new grid ve- 
locity, 
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(15) 



The inner radial flux, F'^^^j^ , is dependent on 
the grid velocity at the inner interface. At the flrst 
radial zone boundary next to the rigid core we im- 
pose both a zero radial velocity and grid velocity. 
Thus, equation (jl5p can be solved recursively from 
the model interior boundary to the surface to de- 
termine the grid velocity at all interfaces. 
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3. COMPUTATIONAL SETUP 

3.1. Starting Model 

The initial model for our adiabatic simulations 
is generated by requiring that it be in hydrostatic 
equilibrium. When this constraint is applied to 
the conservation equations the only terms that re- 
main are the pressure and gravity terms in the 
radial momentum conservation equation. In par- 
ticular, there are no terms left in the internal en- 
ergy conservation equation, and thus no equation 
to solve for the energy structure. To provide this 
information, an energy profile was gener ated from 
another stellar modeling code ROTORC ( Deupre j 
[1990,) 

and energies were interpolated in log(Mr) 
to cell centers. Once we impose the energy dis- 
tribution, we can simultaneously solve the radial 
hydrostatic equilibrium finite difference equation 
and the equation of state for the pressure and 
density structure of the model given the spacing 
of the independent variable Mr- The radius is 
determined from the volume required to produce 
the calculated density from the mass of the shell. 
No convective model is included in the starting 
model because RR Lyrae do not have extensive 
convective regions to affect the structure. To in- 
duce pulsation a radial velocity profile from the 
linear, non-ad iabatic, radial pulsation code, LNA, 
(|Castojll97l"h . modified to allow a gamma law gas, 
is imposed so that the model pulsates around the 
equilibrium point in either the fundamental or the 
first overtone modes. 

3.2. The Grid and Numerical Details 

The simulation volume is separated into cells 
bounded by intersecting surfaces. These surfaces 
are defined at constant values of the three indepen- 
dent variables Mr, 9, and 0. Dependent quantities 
p, E, and P are defined at cell centers and depen- 
dent quantities r, Vr, vor, vg, and are defined 
at appropriate cell interfaces. The models used 
for testing have 107 radial, 10 theta, and 10 phi 
zones. The inner 10 radial zones are handled in 
ID as discussed in ^ The zone number at which 
the switch between ID and 3D is made is chosen 
by the user. For the test cases used in this pa- 
per the total mass of the star was 0.575 Mq, with 
an initial mass spacing of 4.5 x 10~^ Mq at the 
surface, and increasing by 10% each shell into the 
star. Both the 6 and (p zones have a spacing of 



1°, so that the total simulation volume covers 100 
square degrees. 

The equations outlined in H2.1\ are in differen- 
tial form and are approximated by appropriate 
finite difference expressions. Spatial differentials 
are approximated by differences between quanti- 
ties at either cell centers or cell interfaces depend- 
ing on whether the quantity being updated in time 
is interface centered or cell centered, respectively. 
Temporal differentials are approximated by dif- 
ferences between the current grid state and the 
updated grid state divided by the time step. At, 
computed as a fraction of the minimum time step 
allowed by the Courant condition for the model as 
a whole. This then allows us to explicitly solve for 
the updated grid state given the current grid state 
and the time step. 

Equation ()10p is written in finite volume form 
with fluxes defined at cell faces. The veloci- 
ties required for these fluxes are already inter- 
face centered; however, the densities are not. 
In general, quantities that are needed at inter- 
faces but defined at cell centers, and quantities 
that are needed at cell centers but defined at 
interfaces are approximated by averages of ad- 
jacent quantiti es. We have us ed artifi cial viscos- 
ity given bv (Ivon Neumann fc Richtmvei 1950l : 
Richtmver &: Mortonlll967[ ) to smooth out shocks 
with a threshold velocity of one-hundredth of the 
local sound speed for turning on the artificial vis- 
cosity and have used weighted donor cell to stabi- 
lize advection terms, with a weight of 0.1 on the 
upwind terms and 0.9 for centered terms. 

3.3. Order of Calculation 



The order of calculation follows that of lDeupree 
(|l977al ) with a few minor modifications. We start 
by updating the three velocities using equations 
(m m and[5l) from time n — 1/2 to n -I- 1/2 using 
quantities at n {p, (p), r, and P) and quantities 
at n — 1/2 {vr, Vro,vg, and v^). Next the grid ve- 
locity is calculated at time n -f 1/2 using equation 
()15|) working from inner boundary of the model to 
the surface in a recursive manner. The updated 
radius is computed with equation ([7|). The den- 
sity is updated from n to n -I- 1 using the equation 
for mass conservation (eq. |10j). with quantities 
at n [p and r), and quantities at n -I- 1/2 [vr, Vro, 
vg, and v^). The energy is updated in a similar 
manner. The equation of state then allows us to 
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compute the pressure at the new time step from 
the updated density and specific internal energy. 

3.4. Parallelism 

The code we have developed to perform these 
calculations has been named SPHERLS (Stellar 
Pulsation with a Horizontal Eulerian Radial La- 
grangian Scheme). SPHERLS has been designed 
from the beginning to allow for parallel calcula- 
tions using MPI protocols. The parallel design 
allows for domain decomposition in all three di- 
rections with the ability to vary the number of 
ghost cells (used to express the boundary condi- 
tions of the local domain) copied from other pro- 
cessors. Note that boundary conditions in this 
sense are not the global boundary conditions of 
the calculation but only the information required 
from other processors to be able to perform the 
calculations on the processor in question. Equa- 
tions (O, O, ©, ©, (Uni) depend on only local 
quantities and are easily applied to the local grids 
on each processor. The equation to calculate the 
grid velocity (eq. [15| ) requires information across 
all j and k space. Using this equation with domain 
decomposition in the j and k directions would re- 
quire additional message passing which has not yet 
been implemented, and thus currently limits do- 
main decomposition to the radial direction only. 
This could change in the future and may become 
helpful for optimizing calculations for larger hori- 
zontal grids. 

During program initialization, each processor 
can be assigned different equations to solve on 
their local grid, which allows one to divide the 
computational domain up into a 3D and a ID re- 
gion. The ID region is composed of only a single 
zone at each i spanning all of j and k space. Quan- 
tities are volume averaged from the 3D region to 
the ID region to be used as boundary conditions 
for the ID region; while values in the ID region 
are copied across j and k space to the 3D region 
to be used as boundary conditions there. 

To date, only exploratory time trials have been 
performed because future additions to the code 
(including both an implicit solution to the energy 
equation with radiation diffusion and an eddy vis- 
cosity sub-grid-scale model) will impact the tim- 
ing results significantly. At present, a calculation 
with 97 X 10 X 10 3D zones and 10 ID zones for 1 
million time steps (10 million seconds, or approx- 
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Fig. 2. — This figure shows a two-dimensional slice 
at constant (3'^^ (j> zone) for Callll. The color 
scale is (pi,j,k — {p)i)/{p)i) and vectors show the 
difference between the radial velocity and the grid 
velocity added vectorially with the 9 velocity. The 
slice is at 7225 s into the calculation. Cells exterior 
to the white lines on the left and right sides are 
used to express the horizontal periodic boundary 
conditions. This figure shows only the outer 30% 
of the stellar radius, while the total simulation is 
of more than 85%. 

imately 178 fundamental mode periods) takes 2 
processors 12''22™; 4 processors 5^*25"'; 8 proces- 
sors 3''28™; and 16 processors 2^*29'". At larger 
numbers of processors with the current gridding, 
the overhead from MPI begins to negate any addi- 
tional benefits. Thus, for this gridding 16 proces- 
sors represents the "sweet spot". At larger hor- 
izontal grid sizes the "sweet spot" will likely be 
pushed to larger numbers of processors. 

4. TEST CASE RESULTS 

Five adiabatic test calculations have been per- 
formed using the method outlined above. The 
first calculation (Call) was of a static stellar model 
with all velocities set to zero, which was used to 
test that the starting model was indeed gener- 
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Fig. 3. — This figure shows the geometry of the 
torus used to define the toroidal velocity pertur- 
bations, xi and X2 are two radii used to define the 
equation of torus, a and f5 are two angles used to 
define a point on the surface of the torus. The 
upper panel shows a top down view of the torus, 
while the lower panel shows the side view along 
the slice axis indicated in the top panel. The dis- 
tance from an arbitrary point P to the surface of 
the torus is given by d. 

ated in equilibrium consistently with respect to 
the hydrodynamic finite difference equations, and 
that SPHERLS's finite difference and finite vol- 
ume representations of the hydrodynamic equa- 
tions are hydodynamically stable. The second cal- 
culation (C alll) was a spherical blast wave (e.g. 
Sedo3[l95i), used to test that the code could han- 
dle strong shocks and check it against an analytical 
solution. The third calculation (Callll) was of a 
low amplitude radial pulsation (1 km s^^ initial 
surface velocity) in the fundamental mode with a 
horizontal velocity perturbation to break spheri- 
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Fig. 4. — This figure shows a two-dimensional slice 
at constant cj) (6*^ (j) zone) with vectors showing 
the difference between the radial velocity and the 
grid velocity added vectorially with the 9 velocity. 
This plot is from CalV at 9031 s into the calcu- 
lation. At later times the initial toroidal velocity 
perturbation has spread through out the model, 
making its initial form indiscernible. 

cal symmetry. The forth calculation (CallV) was 
of an even lower amplitude radial pulsation (0.1 
km s^^ initial surface velocity) in the first over- 
tone mode again with a horizontal velocity per- 
turbation. These two low amplitude calculations 
were used to compare the calculated periods to 
the linear adiabatic periods, and to insure that 
the scheme worked as expected at low velocities. 
The lower velocity is needed for comparison with 
the linear LNA code which assumes small pertur- 
bations from the static model to calculate periods. 
The fifth calculation (CalV) started with the same 
stellar model as Callll, but with a higher ampli- 
tude pulsation (10 km s~^ initial surface veloc- 
ity) and a toroidal velocity perturbation instead of 
the horizontal velocity perturbation. This calcula- 
tion was used to show that SPHERLS reproduced 
quantities well from one period to the next over 
many periods in the presence of weak shocks and 
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Fig. 5. — This figure sliows a t"wo-dim.ensional slice 
at constant (6"^ zone) with vectors showing the 
difference between the radial velocity and the grid 
velocity added vectorially with the (f> velocity. For 
the same calculation and time as figure 31 



large scale structured motions. 

The static calculation (Call) started with all 
velocities (r, 0, and cj)) set to zero and was com- 
puted for over 20 million seconds (356 fundamen- 
tal periods). The radial velocities in the surface 
zone reached the largest amplitudes. In this zone 
the radial velocity amplitude initially grew from 
zero to a temporal mean of about 3.7 x 10~^ cm 
s~^ within the first 80 periods. This mean was 
maintained for the rest of the calculation with a 
standard deviation of 2.7 x 10~^ cm s~^. The hori- 
zontal velocities remain zero throughout the calcu- 
lation. This is understandable because all terms 
in the horizontal momentum equations are zero 
and remain that way; while the radial momentum 
equation is the balance between two nonzero terms 
and thus subject to round off error. 

To assure that the method is behaving as de- 
signed we checked the mass calculated from the 
averaged density, (p), and the shell volume (which 
is dependent on the radial grid velocity) with the 
mass set as the independent variable. The largest 
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Fig. 6. — This figure shows the radial grid velocity 
for three well separated periods, at six radial shells 
of the model for CallV. 
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Fig. 7. — This figure is similar to figure [6] but for 
(p). The density is well reproduced across a large 
span of periods through the model, indicating that 
there is no drift of the radial coordinate system. 
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relative difference between the calculated mass 
and the independent mass variable for all calcu- 
lations is 4 X 10~^^. This is only two significant 
digits above machine round off, and there are no 
signs of a trend with time. 

To test how well SPHERES handles strong 
shocks and to compare the computed results with 
an analytic calculation we performed a Sedov blast 
wave calculation (Calll). This calculation had 400 
radial zones with an initial spacing of 10 cm and 
10 6* and (j) zones with spacings of 1°. The 400 ra- 
dial zones produces a 40 m radius spherical volume 
for the shock to expand into, and was chosen to 
allow enough volume for the shock to expand into 
over 10 ms, at which time the analytic solution has 
reached a shock radius of 32.7 m. The inner 10 ra- 
dial zones were (as in the adiabatic stellar models) 
treated in ID. The blast was accomplished by set- 
ting the initial energy in the inner 30 zones (i.e., 
a 3 m radius sphere) to 4.18 x 10^^ ergs with all 
other zones having an energy of 1 x 10^ ergs. The 
density was set to 2 g cm^'^ through out the start- 
ing model. A gamma-law gas was used for the 
equation of state, with a 7 of 1.6. All the initial 
velocities were set to zero and the blast was fol- 
lowed for 10 ms. The calculation was compared to 
an analytical solution with a point-source energy 
producing the blast, evaluated at times from 0.5 
ms to 10 ms in 0.5 ms intervals. When comparing 
the extended-source to the point-source solution 
one would expect that at early times (when the 
shock is closer to the source), and at later times 
closer to the initial location of the source, the dis- 
crepancies between the computed and analytical 
solutions should be larger. This is because the 
differences between a non-point source calculation 
and the point source analytic solution will dimin- 
ish as the disturbance moves outward. The blast 
radii computed by SPHERES matched those form 
the analytic solution to within 7.5 cm at all times. 
The best match of shock radii (within 1.7 cm) oc- 
curred later in the calculation at a time of 10 ms, 
while the worst match (7.3 cm ) occurred much 
earlier in the calculation at 3.5 ms. The computed 
velocity, density, and pressure profiles were also 
compared to the analytic solution, but because of 
the extended source, only the zones outside the 
initial explosion source were compared. The root 
mean square of the fractional error in velocity, den- 
sity, and pressure was less than 3%, 8%, and 5% 



respectively in the last half of the calculation (5 ms 
to 10 ms). In the first half of the calculation (0.5 
ms to 5 ms) the root mean square of the fractional 
errors are a bit larger, mostly due to the difference 
between using an extended source in the calcula- 
tion versus a point source in the analytic solution 
and are within 8%, 15%, and 11% for the velocity, 
density, and pressure respectively. The radial pro- 
files of the velocity, density, and pressure fit quite 
well without any outlying points. 

Because the calculations are adiabatic, we ex- 
pect the pulsation to neither grow nor decay and to 
be reproducible from one period to the next. This 
should provide a good test to verify that our nu- 
merical algorithm functions as desired over many 
periods. Both the low amplitude fundamental and 
first overtone pulsation (Callll and CallV respec- 
tively) had a horizontal velocity perturbation im- 
posed on them to break spherical symmetry, by 
setting specific values of vg and at a central 
horizontal zone located at 90% of the total radius 
(18 zones in from the surface of the 107 radial zone 
models). The velocities were directed horizontally 
out of the zone through sides ^j±i/2 and Ai^±i/2 
(see figure [T]). The magnitude of these horizontal 
velocities was taken to be half of the initial radial 
velocity at this radial location (0.3 km s~^ and 
0.03 km s~^ for Callll and CallV respectively). 
Eigure[2]shows a two-dimensional slice at constant 
(j) of Callll slightly after the initial conditions. The 
slice is at 7225 s into the calculation (relatively 
early in the 1 x 10^ s calculation) and shows the 
disturbance resulting from the horizontal velocity 
perturbation as well as it's location and geometry 
with respect to the rest of the model. 

The period of the fundamental mode is 56178 
s for the low amplitude calculation (Callll), and 
compares well with the period calculated from 
ENA of 56114 s. There is less than 0.12% dif- 
ference between the periods of the two codes. The 
first overtone model (CallV) was found to have 
a period of 38911 s and compares with the ENA 
period of 39522 s, producing less than a 1.6% dif- 
ference. 

In addition to the horizontal velocity perturba- 
tion we explored in Callll and CallV we also ex- 
plored a velocity perturbation that is more struc- 
tured over a larger scale (CalV). To create this 
model we started with the same structural model 
and radial velocity profile as Cal HI, this time how- 
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ever, using a surface amplitude of 10 km s~^. On 
top of the radial velocity profile we added a veloc- 
ity perturbation in the shape of a torus (see figure 
[3] for torus geometry). The velocity perturbations 
were taken to be constant on the surface of the 
torus (the two circles in the lower half of figure 
[3|) and parallel to the surface of the torus. By lo- 
cating the closest point (defined by angles a and 
/3) on the surface of the torus to the point P, the 
distance d in figure |3] can be calculated. Then a 
Gaussian centered on the surface of the torus with 
a maximum amplitude of 5 km s""'^ is evaluated at 
d providing the velocity magnitude. The FWHM 
of the Gaussian is chosen so that the velocity per- 
turbations do not overlap the other parts of the 
torus, and so that the velocity perturbations are 
still reasonably strong a zone or two away from 
the surface of the torus. The direction of the ve- 
locity is taken to be parallel to the surface of the 
torus at the locaiton closest to P. The velocity 
magnitude is then broken into r, 9, and (jj compo- 
nents. The result of applying this perturbation is 
shown in figures S] and O These figures show slices 
through the center of the model at constant 9 and 
(p respectively, 9031 s into the calculation and in- 
dicate that 9 and (f> directions behave identically. 

The high surface velocity model (CalV) results 
are presented in figures [6] and [T] Figure [6] shows 
that the radial grid velocity throughout the model 
is well reproduced at periods 50, 110, and 170. 
Figure [7] shows that (p) is well reproduced over 
a large number of periods and d oes no t show 
the drifting apparent in (Dcuprcci Il977a^ . The 
fact that the calculations of dependent quantities 
are reproduced over many periods shows that the 
scheme is working as desired and we expect to cal- 
culate full amplitude solutions of pulsating vari- 
able stars in the future. 

A low radial surface velocity model(l km s~^ 
surface velocity) that was spherically symmetric 
and did not use artificial viscosity has a very con- 
stant peak kinetic energies per period, with out 
any long term detectable growth or decay rates 
larger than 1 x 10""% per fundamental period. 
While Calm, that had a horizontal velocity per- 
turbation, had a growth rate of 4.3 x 10^"*% per 
fundamental period. In the high radial surface ve- 
locity calculation (CalV) artificial viscosity is re- 
quired in the momentum equations (eq. [31 HI and 
[5]) to reduce the very steep pressure gradients at 



shocks. If it is omitted when pulsation amplitudes 
exceed the local sound speed it ultimately leads 
to negative densities and energies. One might ar- 
gue that including the artificial viscosity in the en- 
ergy equation would produce a non-adiabatic cal- 
culation, since the inclusion of the artificial viscos- 
ity raises the internal energy more than otherwise 
when the volume is decreasing. In CalV if the ar- 
tificial viscosity is included in the energy equation 
we find that the peak kinetic energy decays at a 
rate 0.36% per fundamental period. If it isn't in- 
cluded the peak kinetic energy decays at a rate of 
0.13%. These decay rates are dependent on the 
inclusion or omission of artificial viscosity in the 
energy equation, and may affect the amplitude of 
full amplitude solutions to some degree. This is 
not merely a concern for the present code, but for 
all non-linear hydrodynamics codes that use arti- 
ficial viscosity. 

5. CONCLUSIONS AND NEXT STEPS 

We have developed a numerical algorithm and 
the computer code, SPHERLS, which is able to 
follow 3D adiabatic radial pulsations for many pe- 
riods when spherical symmetry is broken. This is a 
necessary step for following the convectivc motions 
and radial pulsations of stars. We have shown that 
SPHERLS maintains hydrostatic equilibrium to a 
high degree over 356 fundamental periods. The 
radial Lagrangian algorithm for maintaining con- 
stant mass in the radial shells is effective over 178 
periods to within one or two digits of machine 
round off, with no signs of any particular trend 
either decreasing or increasing. 

We have found that SPHERES reproduces 
both the fundamental and first overtone modes 
of the linear adiabatic code LNA reasonably well 
(<0.12% for fundamental and < 1.6% for first over 
tone). The velocities as well as other dependent 
quantities (e.g. the horizontally averaged density) 
are reproducible over many periods when spherical 
symmetry is broken, again indicating that our ra- 
dial Lagrangian scheme is performing as designed. 

Adiabaticity is maintained extremely well at 
low radial velocities when artificial viscosity is not 
included. However at higher radial velocities some 
kinetic energy is converted into internal energy 
via the artificial viscosity required to smooth out 
shocks. This should be kept in mind while using 
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any non-linear hydrodynamics code employing ar- 
tificial viscosity to compare full amplitude solu- 
tions with observations. 

In order to calculate full amplitude pulsating 
models to compare with observations a few addi- 
tions must be made to SPHERLS. A more realis- 
tic equation of state and radiative Rosseland mean 
opacities must be included and radiation diffusion 
must be added to the energy equation. The lat- 
ter is expected to require an implicit integration 
of the energy equation, at least near the surface, 
because optically thin zones would require a very 
small time step based on the speed of light and not 
the speed of sound. We expect to use the current 
explicit solution to the energy equation deeper in 
the envelope where the mean free path of photons 
is much less than a computational zone, keeping 
the calculation time down. Finally, we will add 
the subgrid scale terms required in a large eddy 
simulation for treating turbulent convection in the 
ionization zone. 

These calculations were performed with ACEnet 
computational resources. ACEnet, a part of Com- 
pute Canada, provides academic high performance 
computing to Atlantic Canada. CMC is supported 
in part by an NSERC Discovery Grant to RGD 
and in part by an ACEnet fellowship. 
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